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ABSTRACT 

We present Karl G. Jansky Very Large Array (VLA) observations of the CO(2-l) emission in the 
z = 4.05 submillimeter galaxy (SMG) GN20. These high-resolution data allow us to image the 
molecular gas at 1.3 kpc resolution just 1.6 Gyr after the Big Bang. The data reveal a clumpy, 
extended gas reservoir, 14 ± 4 kpc in diameter, in unprecedented detail. A dynamical analysis shows 
that the data are consistent with a rotating disk of total dynamical mass 5.4 ± 2.4 x 1O 11 M . We 
use this dynamical mass estimate to constrain the CO-to-H2 mass conversion factor (aco), finding 
qco = 1-1 ± 0.6 M Q (K km s _1 pc 2 ) -1 . We identify five distinct molecular gas clumps in the disk 
of GN20 with masses a few percent of the total gas mass, brightness temperatures of 16-31K, and 
surface densities of > 3, 200— 4, 500 x (aco/0.8) M pc -2 . Virial mass estimates indicate they could 
be self-gravitating, and we constrain their CO-to-Bb mass conversion factor to be <0.2-0.7 Mq(K 
km s _1 pc 2 ) -1 . A multiwavelength comparison demonstrates that the molecular gas is concentrated 
in a region of the galaxy that is heavily obscured in the rest-frame UV/optical. We investigate the 
spatially-resolved gas excitation and find that the CO(6-5)/CO(2-l) ratio is constant with radius, 
consistent with star formation occuring over a large portion of the disk. We discuss the implications 
of our results in the context of different fueling scenarios for SMGs. 
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galaxies: star formation 



1. INTRODUCTION 



Sub millimeter-luminous galaxies (SMGs; B lain et al.l 
120021 ) are dusty, gas-rich, high-z galaxies that were re- 
vealed in the first extragalactic surveys using SCUBA 
and MAMBO (e.g ISmail et al 



997t|Hughes et al 
1999; iBlain et al.l 



1999; 



1998- ^ 



Bargcr et al. 1998; Ea les et ah. 

Bertoldi et al.ll200aiGreve et al.ll200i . Their huge far 
infrared (FIR) luminosities (^10 13 L©) are believed to 
be primarily driven by intense (^10 3 M^, yr -1 ) star for- 
mation (e.g., lAlexander etall 120051 Mim , adding to 
their already significant st ellar masses (e.e.. iBorvs et all 
120051: lHainline et~aTll201lD . Together with FIR luminous 
quasars, they are generally thought to evolve into the 
giant ellipticals we see in the local universe. They are 
therefore critical for understanding early-type massive 
galaxy formation. 

SMGs are a relatively rare phenomenon, with typ- 
ical space densities of 10~ 5 -10~ 6 Mpc -3 . Their rar- 
ity may be partly because their e nhanced SFRs ar e 
necessarily short-lived (< 100 Mvr: iGreve et al.l [20051. 
or else the galaxies would grow too large, too fast. 
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The highest volume density of radio-selected SMGs 
occurs at z ~ 2 — 3, indicating that they peak si- 
multaneously with the peak epoch of star formation 
(|Chapman et al.ll2003af) . Howeve r, a number of recently 
discovered higher-r ed shift SMGs dSchinnerer et al.ll2008l: 
D addi et al.ll2009bjlal: IRiechers et al.l l2010t iCoppin et al. | 
| 2010tlCapak et al.lboilHWardlow et al.ll2011b ICox et al.l 
20111 ) as well as several studies based on statistical ar- 



guments (|Greve et al.l 120081: iPenner et aT]|201lL and ref- 
erences therein) may indicate that a high-redshift tail 
does exist and is able to account for the population of 
old, massive elliptic als already in place at z ~ 2 — 3 
(jDaddi et al.H2009bD . 

Many SMGs are believed to be starburst-dominated 
major mergers (e.g.. iChapman et al.ll2003bt .lEngel et al.l 
l2lnalNarayanan et al.ll2010HHavward et al.ll20111l2012f) . 
This would make SMGs the high-redshift analogs of 
Ultra-Luminous Infrared Galaxies (ULIRGS) in the lo- 
cal universe. Indeed, there is direct evidence for multi- 
ple CO components and/or disturbed kine matics in some 
SMGs, supporting the merger picture (e.g . .lTacconi et all 
2008HEngel et al.ll20ldHlvison et aTJl2011t IRiechers et al.l 
2011bD . 



Recently, it has been suggested that other mechanisms 
that drive extreme star formation rates may also be at 
play. In particular, a scenario has been put forward 
in which the star formation in massive, high-redshift 
galaxies is driven by cold mode acc retion (CMA; e.g. 
iKeres et al.l 120051: iDekel et al.l l2009aHH ). CMA-driven 
galaxies are constantly forming stars at high rates, and 
the star formation is sustained by smooth infall and ac- 
cretion of gas-rich material. This process can result in 
elliptical galaxies because the streams that feed star for- 
mation can cause the disk to break up into giant clumps 
if they have a high enough gas fraction and degree of 
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turbulence. The clumps th en potentially migra te inward 
and merge into a spheroid (|Dekel et al.ll2009aHb1 ). 

The CMA phenome non has been extended to SMGs by 
a number of authors (Fardal et al.l 120011: iFinlator et al.1 



120061 iDave et al.ll2010l ICarilli et alJ 12010( 1 . In this sce- 
nario, SMGs are massive galaxies sitting at the centers of 
deep potential wells and fed by smooth accretion. They 
can be thought of as super-sized versions of normal star 
forming galaxies, providing an alternate formation sce- 
nario to the merger- induced model. By identifying sim- 
ulated SMGs as the most rapidly star-forming systems 
that match the number densities of SMGs, this theory 
has been successful at explaining some key SMG proper- 
ties, including the ir stellar masses and clustering scales 
(|Dave et al.ll20Toh. 

According to iShapiro et all ((2008), the best way to 
distinguish between CMA and a gas-rich merger model 
is through observations of the gas dynamics and dis- 
tribution. A well defined disk with a smooth rota- 
tion curve would be indicative of CMA, while tidally 
disturbed gas, with a very high Tb starburst nucleus, 
would point towards a merger. This is simply due to 
the redistribution of angular momentum. The large- 
scale gravitational torques induced by gas-rich major 
mergers are efficient at re moving angular momentum 
ijBarnes fc Hernquistl 11996(1 . thereby funneling the cold 
molecular gas into the galaxy's center and producing a 
nuclear starburst. Indeed, in ULIRGs - the canonical 
merger sce nario - the gas is concentrated in the central 
krjc (e.g., iDownes fc Solomon! 119981 : iBrvant fe; Scovillel 
1999). The gas-rich, star-for ming disks in Arp22 
have a size of only ^100 pc ( Sakamoto et al.l |2008). 
In simulations where SMGs result from mergers, their 
gas reservoirs show slig htly larger extents of a few kpc 
((Narayanan et al.ll2009( l. 

Conversely, the z ~ 2 massive star forming galax- 
ies (SFGs) which have been proposed to be due 
to smooth accretion are large disk galaxies, with 
star f ormation occuring at large radii or situated in 
rings flElmegreen fe Elm eerecn 2006; iGenzel et al.l 120011 : 
lElmegreen et al.1 12009( 1. Rest-frame UV/optical imag- 
ing shows that, unlike low-redshift galaxies, their 
disks tend to be broken into multiple giant c lumps 
of ~1 kpc and 10 9 M w (lElmegreen et al.l 12005 
Elmegreen fc Elmegreenl 20051 : iForster Schreiber et al.l 
20061: IGenzel et all 120081 1. More recent work has 



also imaged such galaxies in millimeter /CO em ission 
((Tacconi et al.l l20int ISwinbank et "all l2?TToL [201lll . find- 
ing evide nce for clumpy CO e mission extending over the 
disk fe.g. lTacconi et al.|[20Tol l. 

Observations of the morphology and kinematics of 
molecular gas in SMGs may therefore shed light on 
the p hysics behind the in tense star formation. Previ- 
ously, ICarilli et al.1 ((2010( 1 presented VLA observations 
of CO(1-0) and CO(2-l) emission in the z = 4.05 
SMG GN20, the brightest SMG in the GOODS-N field 
(|Pope et al.lf2006h and one of three SMGs in what ap- 
pears to be a massive z ~ 4 proto-cluster ((Daddi et al.1 
2009b). These low- J transitions are of particular inter- 
est because they trace the cold molecular gas thought 
to make up the bulk of the molecular gas in these sys- 
tems. Although their high-resolution CO(2-l) imaging 
showed evidence that the gas in GN20 was well-resolved, 
their observations suffered from some severe spectral lim- 
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Fig. 1. — CO(2-l) spectrum of GN20, binned into channels of 78 
km s — - 1 . The single Gaussian fit to the data is shown by the dotted 
line, and the double-horn fit is shown by the dot-dashed line. The 
bar at the bottom of the plot shows the velocity range averaged 
over to make Figure [4] The model fit results from the dynamical 
modeling discussed in Section 14.1.21 

itations. In particular, the observations utilized the old 
VLA correlator, with a total bandwidth of only 100 MHz 
(650 km s _1 ), causing the line profiles to be truncated 
on both sides. The observations were also done in con- 
tinuum mode, resulting in no information on kinematics. 

We therefore obtained over 120 hours of time on the 
Karl G. Jansky Very Large Array (VLA; iPerlev et al.1 
1201 If) to image the CO(2-l) emission in the GN20 
field in the B- and D-configurations. The ~20 hours 
of lower-resolution D-array data were presented by 
ICarilli et al.1 ((201 1( 1. Here, we present the full dataset 
(B+D-configurations) on GN20, providing greatly im- 
proved spatial resolution and image fidelity. 

We begin in Section [2] by describing our new VLA ob- 
servations and data reduction of GN20. CO maps and 
the derived gas mass are presented in Section [3] Sec- 
tion [4] constitutes our analysis, including a dynamical 
analysis (|4.1I) , the definition and properties of individual 
gas clumps (|4.2|) . a multiwavelength comparison (|4.3|1 . 
and an analysis of the spatially-resolved gas excitation 
(|4.4|) . We discuss the implications of our results on the 
nature of GN20 in Section [5l and we end with a sum- 
mary of our conclusions in Section [6] Where applicable 
we assume the standard A cosmolog y of Hp = 70 km s -1 
Mpc" 1 , fl A = 0.7, and Q M = 0.3 ((Spergel et al.l [20031 
I2007f) . At a redshift of z = 4.05, 1" corresponds to 7 
kpc. 

2. OBSERVATIONS & DATA REDUCTION 

We observed the CO(2-l) transition toward the GN20 
field as part of VLA key project AC974. The project 
was awarded 96 hours in B-configuration (baselines up 
to 10 km) and 28 hours in D-configuration (baselines 
up to 1 km), for a total of 124 hours. The observations 
were dynamically scheduled and took place in March- 
April 2010 (D-configuration), and February-April 2011 
(B-configuration) . 

The CO(2-l) line (rest frequency v = 230.5424 GHz) 
is redshifted to v = 45.655 GHz at z ~ 4.05, requiring 
the Q-band. The primary beam is ~1' (FWHM) at this 
frequency, and the pointing center was chosen to be 10" 
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Fig. 2. — VLA CO(2-l) in GN20 in 78 km s _1 channels. Increasing channel numbers indicate increasing frequency (decreasing velocity), 
and channel 6 corresponds to the central frequency derived from the spectrum (Figure [TJ. The left panel is at 0.19" resolution, and the 
right panel shows the same channels tapered to 0.38" resolution to increase the S/N ratio. The cross shows the peak of the 1.4 GHz 
counterpart at 1.7" resolution (Morrison ct al. 2010). The rms noise values are 59 fijj beam - and 70 fijy beam , and the contours are 
given in steps of la starting at ±2<r. 



west of GN20 so that GN20, the nearby SMGs (and fel- 
low proto-cluster members) GN20.2a and GN20.2b, and 
the z = 1.5 galaxy BzK-21000 would all fall within the 
70% sens itivity radius of th e primary beam (data pre- 
sented in lHodge et al.l l2012). All images have been cor- 
rected for the response of the VLA primary beam. We 
centered the two 128 MHz IFs at 45.592 GHz and 45.720 
GHz, for a total bandwidth of 246 MHz or 1600 km s" 1 
(taking into account the overlap). Each of the two IFs 
had 64 channels, resulting in an instrumental velocity 
resolution of 13 km s _1 . Observations were taken in full 
polarization mode. 

We used fast swi tching phase calibration 
(|Carilli fc H oldawav 1999), with a cycle time of 5 
minutes. VLA calibrator J1302+5748 served as the 
phase calibrator. The quasar 3C286 was used to deter- 
mine the absolute flux density scale and for bandpass 
calibration. The B-configuration data were reduced 
using the Astronomical Image Processing System 
(AIPS), and the D-configuration data were reduced in 
the Common Astronomy Software Applications (CASA) 
package. During the reduction, we discarded data taken 
during times of poor phase stability. After accounting 
for calibration overheads and flagging, the total time on 
source was approximately 50 hours. The data were then 
combined within AIPS for further analysis. 

We imaged the data using the CLEAN algorithm, and 
we cleaned down to 1.5<r in a 2" x 2" clean box around 
GN20. A spectrum was extracted using an aperture the 



same size as the clean box. The best compromise between 
resolution and sensitivity for this dataset resulted from 
using Briggs weighting with a robust parameter of R = 
1.0. This results in an angular resolution of 0.19", or 
1.3 kpc at the redshift of GN20. We will refer to this 
resolution as our 'native' resolution for the remainder of 
the paper. 

3. RESULTS 
3.1. CO (2-1) Detection and Channel Maps 

The CO(2-l) spectrum for GN20 is shown in Figured] 
We used a Gaussian fit to the spectrum at 12 MHz (78 
km s _1 ) resolution to derive a peak flux density of 1.35 
± 0.20 mJy and a FWHM of 730 ± 140 km s"\ con- 
sistent with ICarilli et al.l ([201 ID . The integrated flux is 
1.0 ± 0.25 Jy km s _1 (see Table Q]). The spectrum indi- 
cates a redshift for the galaxy of z = 4.0548 (topocentric), 
consistent with earlier measure ments, including high- J 
transitions (jDaddi et al.l l2009bf) . We detect no contin- 
uum emission from GN20 at our native resolution (0.19" 
with a 3cr limit of 39.5 /iJy and at 1.75" resolution with 
a 3(7 limit of 65.4 //Jy. 

The spectrum shows the clear indication of a double- 
peaked, or at least flattened, structure, the signature of 
which is also evident in a position-velocity diagram cre- 
ated at a later point in the analysis ( Section I4.1.2[) . If 
we instead fit the spectrum with a combination of two 
Gaussians (Figure [TJ dot-dashed line), we find that the 



4 



Hodge et al. 



TABLE 1 

GN20 Observed and Derived Parameters 



Parameter 



Value 



Position (J2000) a 



5 CO(2-l) 



FWHM. 



CO(2-l) 



J CO(2-l) 

T I 

^CO(2-l) 

M(H 2 ) 
M dyn 



12 h 37 m 11.89 s +62°22'11.8" 
4.0548 ± 0.0008 
1.6 ± 0.3 mjy (peak 1) 
1.3 ± 0.3 mjy (peak 2) 
300 ± 90 km s" 1 (peak 1) 
320 ± 130 km s" 1 (peak 2) 
1.0 ± 0.3 Jy km s" 1 
1.6 ± 0.5 x 10 11 K km s" 1 pc 2 
1.3 ± 0.4 x 10 11 x (aco/0.8) M Q 
5.4 ± 2.4 x 10 11 Mo 



* From the 1.4 GHz observations of Morrison et al. 

(|20TO| ) at 1.7" resolution. All other parameters are 

from the study presented here, 
two components have peak flux densities of 1.6 ± 0.3 mjy 
and 1.3 ± 0.3 mjy, and FWHM values of 300 ± 90 km 
s _1 and 320 ± 130 km s _1 , respectively. The integrated 
flux of the combined components is 1.0 ± 0.3 Jy km s -1 , 
consistent with the result from the single-Gaussian fit 
above. We will use the two-Gaussian fit to determine 
the CO luminosity and gas mass in the remaining anal- 
ysis. 

Channel maps of 78 km s _1 width are shown in Fig- 
ure [2] The left panel shows the emission at our native 
resolution (0.19"), and the right panel shows the same 
channels tapered to 0.38" resolution in order to increase 
the S/N ratio. The two outer channels in each case sam- 
ple the continuum on either side of the line. The emission 
first appears in the South of the plotted field, and it shifts 
to the North with increasing frequency. We will analyze 
the kinematics of the system in Section 14.11 below. 

3.2. Molecular Gas Mass 
From the two-Gaussian fit to the spectrum of GN20, 
we derive a CO luminosity of ^^0(2-1) = ^ ^ 0.5 
x 10 11 K km s _1 pc 2 . This implies a molecular gas 
mass of M(H 2 ) = 1.3 ± 0.4 x 10 11 x (a co /0.8) 
M Q (see Table 1) assuming the s tanda rd relation- 
ships fr om ISolomon fc Vanden Bout] ()2005D . As justi- 
fied by ICarilli et al.l (120101) . we assumed thermal exci- 
tation for the extrapolation from CO(2-l) to CO(l- 
0), and we used a CO-to-H 2 conversion factor of aco 
0.8 M Q (K km s^ 1 pc 2 )" 1 , the value typically as- 



sumed for ULIRGs and SMGs dDownes k, Solomon! 



Solomon & Vand en Bo ut 2005| iTacconi et al 
2008|). In Section|4T3]below, we use our dynamical mass 



1998 



2006 



estimate to put constraints on aco for this system. 

4. ANALYSIS 

4.1. Dynamical Analysis 

4.1.1. Moment Maps 

Given the high-quality of the data, we created moment 
maps by following the typical approach us ed for HI and 
CO observations of nea rby galaxies (e.g. IWalter et al.l 
120081: ILerov et al.l 120091) . This approach is useful for 
high resolution datasets such as this (0.19"), where real 
emission can be resolved out in a moment map with a 
single, global S/N cut. To address this fact, this tech- 
nique involves using the data, tapered to a lower res- 
olution (hence recovering more extended emission), as 
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Fig. 3. — CO(2-l) 0th (top) and 1st (bottom) moment maps 
for GN20 at a resolution of 0.19". The 0th moment map (i.e. 
integrated intensity) has a peak S/N of 6, and the contours shown 
start at (and are in steps of) 15.5 mjy km s — 1 . Contours for the 
1st moment map (i.e., intensity-weighted mean velocity) are shown 
for steps of 100 km s , with the contour between the green and 
orange bands representing the systemic velocity. 

an additional input. For this purpose, we used data 
tapered to 0.38" resolution. This lower-resolution dat- 
acube is used to create a mask (on a per-channel basis) 
of the significant emission, including more diffuse emis- 
sion. These masks are then used to blank the higher- 
resolution data, thereby retaining real, diffuse emission 
in the high-resolution data that would have otherwise 
been blanked. This standard procedure is the best pos- 
sible method to recover diffuse, low S/N emission while 
retaining the highest possible spatial resolution. 

We used the method described above to create the 0th 
moment (integrated intensity) map, shown in Figure [3] 
(top panel). As the masking process is done on a per- 
channel basis, the noise in any given pixel in the map is 
given by y/N x cr c han, where N is the number of channels 
that were integrated for that particular pixel, and c c han 
is the rms noise per channel. We used this information 
to apply an additional S/N cut to the 1st moment map 



A clumpy, rotating gas disk in SMG GN20 



5 



62 22 14.0 

13.5 - 

: 13.0 

\ 12.5 

I 

12.0 
! 11.5 

I 

11.0 - 
10.5 



10. 




O 



VS 3712.1 




a. 



12.0 11.9 11.8 11.7 



12 3712.1 



12.0 11.9 11.8 11.7 
RIGHT ASCENSION (J2000) 



12 3712.1 12.0 11.9 11.8 11.7 



ing; while SMGs typically have relatively compact dis- 
trib utions in the higher - J CO lines (HWHM=2-4 kpc; 
e.g. iTacconi et al.l I2006T) . recent observations of SMGs 
in CO (1-0) show more extended gas res ervoirs (e.g. 
Ivison et alj 120111: IRiechers et all I2011allbl ; llvison et al.1 
20101) . 



Fig. 4. — Velocity— averaged B+D-array CO(2-l) images of GN20 over a bandwidth of 780 km s 1 at three different resolutions (from 
left): 0.38", 0.19", and 0.14", as given in the lower left corner of the maps. The cross shows the peak of the 1.4 GHz counterpart and the 
size of the cross indicates the 1.4 GHz beam at 1.7" resolution (Morrison ct al. 2010). The rms noise values are 25.6 fijy beam -1 , 19.0 
/ijy beam -1 , and 25.0 /tjy beam -1 , and the contours are given in steps of la starting at ±2ct. 

(intensity- weighted velocity; Figure [3] lower panel). In 
particular, we required S/N > 3 for the 1st moment map, 
since any remaining unmasked noise will have a large 
effect on the resulting velocity field. Although some noise 
is still present in the outskirts, a clear velocity gradient 
is apparent across the disk. 

As a verification of the blanking process described 
above, we show unblanked velocity-averaged maps of 
GN20 in Figure 2] at three different angular resolutions. 
The left panel has been tapered to 0.38", the middle 
panel is at the native resolution (0.19"), and the right 
panel uses Briggs weighting with R=-0.5 to reach 0.14" 
resolution (1.0 kpc at z = 4.05). In contrast to the mo- 
ment analysis described above, these maps were made 
by simply averaging over 780 km s _1 (the region indi- 
cated in Figure [1]), without first blanking the noise on a 
per-channel basis. The gas distribution appears slightly 
different than that seen in the 0th moment map because 
the channels without emission are included in the aver- 
age. The highest-resolution map confirms that the emis- 
sion is not concentrated in just one or two strong peaks, 
but rather spread out over a larger area. 



The total flux density (averaged over 780 km s -1 ) in 
these velocity-averaged B+D-array maps is 633 ± 67 
/ijy For comparison, the total flux densities in the B- 
array only (extended configuration) and D-array only 
(compact configuration) maps (not shown) are 365 ± 57 
/xJy and 813 ± 65 /xJy, respectively. The total flux den- 
sity in the 0th moment map (800 ± 74 /xJy) is consistent 
with that in the D-array only map, confirming that the 
recovered emission in the 0th moment map is real. The 
blanking process described above therefore allows us to 
achieve the native resolution while still recovering the 
diffuse emission present on larger scales. 

We estimate the size of GN20's gas reservoir from the 
0th moment map (Figure [31 top) . Defining the radius 
as the maximum radial extent of the resolved CO(2-l) 
emission in the map, and conservatively assuming an un- 
certainty in the measurement of ^30%, we derive a ra- 
dius of 1" ± 0.3", equivalent to ~7 kpc (± 2 kpc) at 
z = 4.055. The total diameter of the source in CO(2- 
1) is therefore ~14 ± 4 kpc. The large extent of the 
molecular gas reservoir is not unlike what has been seen 
in some other (lower-z) SMGs in low-excitation imag- 



4.1.2. Dynamical Modeling 

For the dynamical modeling of GN20, we used the 
GALMOD task (part of the GIPSY package). GALMOD 
creates a three-dimensional model using a set of input 
parameters, and it then convolves the model cube to the 
spatial/spectral resolution of the data for comparison. 
We used an input data cube with a spectral resolution 
of 26 km s , and we tapered the data to an angular 
resolution of 0.77" as it was found that higher resolu- 
tions resolved out too much emission to be usable. The 
GALMOD task requires a radial profile as input, which 
(guided by the 0th moment map) we set as an exponen- 
tial radial profile with a slope of -0.4. We used a thin 
disk model and we found that changes in the thickness 
of the disk (within a reasonable range, < few kpc) did 
not result in major changes to the model. 

As the S /N of our data have necessitated using a lower 
resolution (0.77") cube for the modeling, it is not possi- 
ble to measure the intrinsic rotation curve of this source 
directly from the data. Instead, we have assumed a rota- 
tion curve as input and then compared the model to the 
data in position-velocity space. Position-velocity dia- 
grams are shown for three different input rotation curves 
in Figure [5] The columns represent the different rotation 
curves: The left-hand column assumes a flat rotation 
curve, the middle column assumes a rotation curve that 
rises linearly from 0"-0.5" (0-3.5 kpc) to v max , then is 
flat at larger radii, and the right-hand column assumes 
a rotation curve that rises linearly from 0"-1.0" to v max 
(0-7 kpc), then is flat at larger radii. The rows show 
four different slices through the three-dimensional data: 
major axis, taken at a position angle of 25° (top row), 
minor axis (second row), major axis + 30° (third row), 
and major axis + 60° (bottom row). These comparison 
plots demonstrate that, while the resolution of the in- 
put data cube (due to our limited S /N) make it difficult 
to constrain the exact shape of the rotation curve, we 
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Fig. 5. — Position-velocity diagrams for CO(2-l) emission in GN20 for four different slices and three different models (shown as contours). 
The input data cube has a spectral resolution of 0.77" and a spectral resolution of 26 km s _1 . The different slices are shown in the different 
rows and are (from top): major axis (i.e. a position angle of 25°), minor axis, major axis+30°, and major axis+60°. The columns show 
the different models (contours): a flat rotation curve (left), a rotation curve that rises linearly to 0.5" (3.5 kpc) then flattens (middle), and 
a rotation curve that rises linearly out to 1.0" (7 kpc) (right). The velocities on the vertical axis are relative. Greyscalc and thin countours 
show the observed data. The data contours are 55% to 85% in steps of 10%. For all models, the contours are 35% to 80% of the peak 
brightness, in steps of 15%. The left-hand panel shows our best fit model (see text for details). 



can constrain the curve to having reached its fiat part at 
least within 0.5" (3.5 kpc), with a preference for the flat 
rotation curve. 

We therefore find that the velocity field is fully con- 
sistent with a rotating disk with a flat rotation curve 
(as discussed above). By comparing different models to 
the data, we find that the best-fit model is a rotating 
disk with an inclination of i = 30° ± 15°, a maximum 
rotational velocity of v max = 575 ± 100 km s _1 , and a 



dispersion of S = 100 ± 30 km s -1 . Note that the er- 
ror bars are not statistical, but are liberal estimates of 
the uncertainty determined through careful comparison 
between model and data. Note also that deriving the 
dispersion from a spatially and spectrally convolved disk 
model, unlike other mean-weighte d dispersion estima - 
tors, is unbiased by beam smearing (|Davies et al.ll20ill) . 
The spectrum predicted by our best-fit model is a good 
fit to our observed spectrum (Figure [TJ . 
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Fig. 6. — Zoomed-in CO(2-l) 0th (left) and 1st (right) moment maps at 0.19" resolution. The clumps in the final clump sample are 
indicated by the numbered circles. 

with the lowest chi-squared values for their fits resulting 
from more ULIRG-like values. However, their resolution 
did only marginally resolve their sources, and they rely 
on CO(4-3) and even higher-order transitions that may 
be more centrally concentrated (see Section PO)) . 

Our derived dynamical mass allows us to put con- 
straints on the conversion factor for GN20. If we as- 
sume that GN20 is composed of 100% molecular gas, 
then we derive a conversion factor of aco — 3.3 ± 1.8 
M© (K km s" 1 pc 2 ) -1 to account for the total mass in 
the system. This is consistent (within the large uncer- 
tainties) with the Galactic value, but it represents the 
most extreme case. In reality, the total mass will also 
include contributions from stars, dark matter, and dust. 
The stellar component has been estimated to be 2.3 x 
10 11 M Q bv lDaddi etall (|2009b[) based on SED fitting 
t o the ACS throu gh IRAC photometry and assuming 
a IChabrierl ([20031 ) IMF. (Note that the estimate relies 
most heavily on the IRAC data, which is coincident with 
the CO emission.) The contribution from dark matter 
is largely unknown, but observations of high-z galax- 
ies suggest the value ma y be roug hly 25% (jGenzel et al.l 
l2008tlDaddi et al.|[2010h . Using these estimates, and ig- 
noring the contribution from dust, which is thought to 
make up only a tiny fraction of the total mass (2 x 10 9 
: IMaedis et"alll2?jrl . we derive a conversion factoi0 
of aco = 1-1 ± 0.6 M© (K km s" 1 pc 2 )" 1 . This esti- 
mate is in agreemen t with the lim i t aco < 1.0 derived 
recently for GN20 by IMaedis efaLl (pOll using the local 
Mg as /Mdust vs. metallicity relation. 

4.2. Molecular Gas Clumps 
4.2.1. Clump search 

In order to study the properties of the molecular gas 
in GN20 in more detail, we have attempted to identify 
individual gas clumps in the disk of GN20. We used two 
different methods to identify potential clumps, judging 
as reliable only those clumps that were independently 
identified by both algorithms. The first method used the 
AIPS task SERCH. SERCH is a "matched-filter" analysis 

8 Note that the formal uncertainty is larger but includes negative 
values of «co an d i s therefore unphysical. Here, we have assumed 
an uncertainty of ~50%. 



The relatively large value of v max is due to the fairly 
small inclination value; while the quantity v max sin(i) is 
well constrained, the two components are difficult to dis- 
entangle at our angular resolution. Following standard 
procedure for the modeling of low to medium resolution 
HI observations, the final inclination value of 30° was 
chosen to reproduce the resulting ellipticity in the 0th 
moment map. However, we cannot definitively rule out 
larger values of the inclination (within the quoted er- 
ror), and therefore lower values of v max . The uncertainty 
quoted for v max folds in the uncertainty in the inclina- 
tion. 

4.1.3. Dynamical Mass and the GO-to-Hi conversion factor 

Assuming the parameters from the best fit model, we 
derive a dynamical mass for GN20 of 5.4 ± 2.4 x 10 11 
Mq. The uncertainty was estimated assuming la uncer- 
tainties of 100 km s _1 and 2 kpc on the rotational ve- 
locity and radius, respectively. This estimate is based on 
dynamical modeling, making it more robust than previ- 
ous estimates for this source (which also relied on higher - 
J transition lines; IDaddi et al.ll2009bt ICarilli et al.ll2010D . 
We will now use this estimate to set limits on the CO- 
to-H2 mass conversion factor. 

In Section [3~2l we calculated molecular gas masses as- 
suming a CO-to-H 2 conversion factor of aco = 0.8 
M© (K km s" 1 pc 2 )" 1 , the value derived for ULIRGs. 
The mass conversion factor is thought to vary with en- 
vironment, however, depending on several factors in- 
cluding metallicity, excitation , and interstellar medium 
pressure ()Bolatto et al.l 120081 and references therein). 
Generally, it is thought that it m ay decrease for ob- 
jects with large gas surface densities (IDownes fc Solomon! 
11991 [Scmalle et al.l [T997I: ITacconi et al.l 120081) . ranging 
from 0.8 M© (K km s" 1 pc 2 )" 1 for ULIRGs up to ~4.3 
M© (K km s" 1 pc 2 )" 1 for giant molecular clouds (GMCs) 
in the Milky Way. (All values stated here include He- 
lium.) While the ULIRG value of 0.8 M© (K km s" 1 
pc 2 )" 1 is often used for SMGs, there is not yet any 
firm evidence for what the SMG value should be. It 
is poss ible that for t h e mor e luminous SMGs, it is even 
lower. ITacconi et al.l ((2008) find that a Galactic conver- 
sion factor is strongly disfavored for their nine SMGs, 
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Fig. 7. — CO(2-l) spectra of five molecular gas clumps in GN20. The spectra have been binned into channels of 39 km s _1 . The 
dot-dashed line shows the observed coherent velocity at the position of the clump, and the dotted line shows the model coherent velocity. 
The S/N of the spectra is not very high, but still allows us to derive approximate line widths for the individual clumps for the first time. 
The results for individual clumps are summarized in Table [2] 




emplo ying a Gaussian kernel and used in source-finding 
(e.g. iBegum fc Chengalurl 120051 : IKanekar fc Chengalurl 
120051 lAravena et al.l |2012[). and the algor ithm is de- 
scribed in more detail in lUson et all (|199lD . We used 
a cube at our native angular resolution (0.19") and with 
a spectral binning of 39 km s _1 . We then searched for 
emission line sources in the three-dimensional data, opti- 
mizing the search for linewidths between 80-200 km s _1 
in steps of the channel size (^40 km s _1 ). A S/N cut 
of 4.5cr yielded 11 positive peaks and two negative peaks 
in a generous 4"x4" region centered on GN20. Based 
on the number of negative detections, we can expect ~2 
of the positive peaks to be spurious. Ten of the positive 
peaks lie on the disk of GN20, as defined by the Oth mo- 
ment map in Figure [3l and one positive peak (and both 
negative peaks) lie off of the disk. 

To test the robustness of the clumps identified in GN20 
by SERCH, we repeated the search using an independent 
method. The second method is built upon the princi- 
ples of Bayesian inference, and we will briefly describe 
it here. (For a more detailed description, see Lentati 
et al. 2012, in preparation.) This method involves using 
the multimodal nested sampling algorithm MULTINEST 
(jFeroz fc Hobsonl [2008) to efficiently sample the poste- 
rior distribution by fitting a simple parametric model to 
the data cube described both spatially and spectrally by 
a Gaussian. The spatial model parameters are fixed to 
match the clean beam, since the clumps are not expected 
to be resolved. Spectrally, the model is a Gaussian line 
profile described by a profile width <j v and peak ampli- 
tude A p . The priors on a v and A p are set as uniform 
between 20-150 km s _1 and 0.0-0.4 mJy, respectively, 
and are considered independent. 

We then quantified the probability that the fitted 
model describes a real clump by using Bayes ian model 
selection (see e.g. lHobson fc McL achlan 2003). The two 



models considered were: 

(1) - "There is no clump present at that position" 

(2) - "There exists a clump with A p > at that position" 

We evaluated the evidence associated with the pos- 
terior for these competing models, where the evidence 
is the average of the likelihood over the prior. Specifi- 
cally, MULTINEST returns the evidence and associated 
error for each peak in the posterior which can then be 
compared with the evidence for there being no clump 
present. The model corresponding to the posterior mode 
with the greatest probabiliy was then subtracted from 
the data, and the process was repeated until no sources 
above a low threshold probability were found. Individ- 
ual clumps were identified by requiring a minimum dis- 
tance of 0.12" and two frequency channels between their 
centers. If multiple clumps were returned within that 
distance, the clump with the greatest evidence was cho- 
sen. This method returned a list of positive and negative 
clumps, sorted in order of their probability. The nega- 
tive clumps were all very low on the list, with 20 positive 
clumps judged as more probable than the most probable 
negative clump. 

We then cross-matched the lists of clump candidates 
produced by both SERCH and the Bayesian modeling, 
requiring the clump centers returned by the two algo- 
rithms to agree within a beam radius in order to count 
as the same clump. Of the ten potential clumps identi- 
fied on the disk of GN20 by SERCH, six were indepen- 
dently identified by the Bayesian modeling (and these 
were judged six of the seven most probable clumps by 
the Bayesian modeling). The four SERCH clumps that 
were not identified by the Bayesian modeling included 
the two SERCH clumps furthest out on the disk, and 
two SERCH clumps which were closer than the minimum 
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TABLE 2 

Molecular Gas Clumps Properties 



Clump a 


S CO(2-l) 

My] 


FWHM co(2 -i) 
[km s" 1 ] 


S/N 

M 


Tb (rest frame) 
[K] 


M(H 2 ) 
[x 10 9 (aco/0.8) Mq] 


% Total Af(H 2 ) 
[x 10 9 Mq] 


[x 10 9 M ] 


1 


360 ± 100 


90 ± 40 


5.3 


31 ± 8 


4.4 ± 2.2 


3% ± 2% 


<1.0 ± 0.9 


2 


270 ± 60 


150 ± 50 


5.4 


23 ± 5 


6.0 ± 2.4 


5% ± 2% 


<3.4 ± 2.2 


3 


270 ± 60 


140 ± 40 


5.1 


23 ± 6 


5.5 ± 2.1 


4% ± 2% 


<2.9 ± 1.8 


1 


180 ± 50 


160 ± 60 


4.5 


16 ± 4 


4.2 ± 1.9 


3% ± 1% 


<3.6 ± 2.6 


5 


240 ± 70 


140 ± 50 


4.6 


21 ± 6 


4.8 ± 2.3 


4% ± 2% 


<2.7 ± 2.0 



a See Figure [6] 

distance set in the Bayesian modeling. Of the six com- 
mon clump candidates, one was discarded because of an 
unnaturally thin profile (i.e. a single spectral channel). 
The final sample therefore consists of five clumps. These 
are the most probable clumps which were independently 
identified by the two different methods, and therefore 
the five strongest clump candidates. They are labeled 
in Figure H] at their positions on the CO(2-l) 0th (left) 
and 1st (right) moment maps with beam-sized circles 
(1.3 kpc). As these clumps were identified in the three- 
dimensional data, they do not necessarily perfectly align 
with the most significant peaks in the integrated map. 
We also note that, while these clumps were judged to 
be the most reliable clumps, they in no way consistute a 
comprehensive census of the clumps in GN20. 

Spectra for these five clumps are shown in Figure 
where the panel numbers correspond to the clump num- 
bers in Figure [6] Gaussian fits to the spectra are over- 
plotted. The coherent velocity at the position of each 
clump is shown in each panel, and is based both on the 
model velocity field at lower spatial resolution (dotted 
line) and the data velocity field at its native resolution 
(dot-dashed line). In general, the clumps are within 
^100 km s _1 of the (observed) coherent velocity at their 
position on the disk. The fact that they are not cen- 
tered exactly on the observed coherent velocity is due to 
the presence of other mass along the line-of-sight. The 
model velocity field is based on data at a lower spatial 
resolution (0.77") and so does worse (generally) at pre- 
dicting the clump velocities. While these spectra are low 
S /N, they allow us to estimate the properties of individ- 
ual gas clumps in a z=4 galaxy for the first time. 

4.2.2. Clump Properties 

With this data, we attempt to estimate some basic 
properties of the molecular gas clumps. We caution that 
this analysis is based on low S/N data and will need to 
be verified by even higher sensitivity observations. To 
begin with, we use the spectral information to deter- 
mine brightness temperatures of individual clumps with- 
out diluting the values by averaging in velocity. Derived 
brightness temperatures (using the Rayleigh- Jeans ap- 
proximation) range from 3.2-6.2 K and can be read off 
of the right-hand y-axis in Figure [7J After correcting to 
the rest-frame, values for individual clumps range from 
~16Kupto31K (Table [5]). For reference, Planck tem- 
peratures (i.e. the temperatures calculated using the full 
blackbody instead of just the Rayleigh- Jeans approxi- 
mation) are typically ~20% higher. 

Since the clumps appear to be unresolved on the scales 
probed by our observations (1.3 kpc), the brightness 
temperature values are most likely lower limits. Nev- 



ertheless, it is interesting to note that the peak bright- 
ness temperatures derived for the clumps are already 
appro aching (or even exc eeding) the dust temperature 
(33K; Maed is et alJl20ll . If we instead compare to 
the large velocity gradie nt (LVG; e.g. IScoville fc Solomon! 
11974ft gas model fit by iCarilh et al.l 1)20100 . where the 
best-fit model consisted of a lower and a higher exci- 
tation component, the higher excitation component still 
has a kinetic temperature of only 45 K. The change in 
surface area required to make up this temperature dif- 
ference is small. This may indicate that we are close to 
resolving the clumps, though we cannot make this state- 
ment any more quantitative with the current dataset. 

We next fit the linewidths of the clumps. There are 
several factors that contribute to broadening the clump 
linewidths beyond their intrinsic values. One factor is the 
instrumental velocity resolution, which we have removed 
in quadrature from the measured linewidths. Another 
factor is the rotational contribution from the large-scale 
velocity structure of the galaxy. To estimate this con- 
tribution, we have taken the model velocity field and 
measured the velocity gradient across a beam-sized area 
centered on each clump, dividing by two to approximate 
the weighting of the beam. The resulting linewidths are 
listed in Table [5] and have a median value of 140 km 
s _1 . The errors were derived from the Gaussian fits and 
assuming a 30% error on the measured velocity gradient. 

Using the Gaussian fits to the spectra in Figure [71 we 
determined the molecular gas masses of the clumps (Ta- 
ble [2]). For this, we used the same relations and assump- 
tions as in Section 13.21 including an H2-to-CO conver- 
sion factor of aco=0.8 M Q (K km s _1 pc 2 ) -1 . We es- 
timate that the clumps have H2 masses of ^5 x 10 9 x 
(c^co/0.8) Mq, or a few percent of the total gas mass 
each. 

Taking the beam size as an upper limit on the clump 
size, the mass surface densities of the clumps are >3,200- 
4,500 x (aco/0-8) M Q pc~ 2 . Similar surface densities 
have been seen in some compact SMGs with double- 
peake d line profiles (3,000-10,000 M Q pc~ 2 ; lEngel eTaTl 
2010). GMCs in the nearby universe, on the other 
hand, are typically several orders of magnitude (in vol- 
ume) smaller and have surface den si ties approximately 20 
times lower ([Solomon et al.lll987ft . iTacconi et all (|2008f ) 
argue that a major merger is necessary to produce such 
high mass surface densi ties as those seen in SMGs. How- 
ever, |EngeTeFljI] (|2010ft have pointed out that this argu- 
ment alone is insufficient to distinguish a merger. Indeed, 
if we assume that the clumps are spherical and roughly 
the size of the beam (at maximum), then we derive (av- 
erage) volume densities of only 70-110 c m" 3 , consistent 
with the mean density of local GMCs (Sol omon et al.l 
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[1981 . The dense cores of GMCs show much higher val- 
■T0 4 cm- 3 (e.g. iGoldsmith et al.lll997f ). and the 



ucs, 



same may also be true for the cores of the clumps in 
GN20. The mean H2 volume density of GN20 as a whole 
is only ~15 cm -3 . 

We estimated the virial mas ses of the clump s using the 
isotropic virial estimator (e.g. lErb~e t al. 200y): 



G 



(1) 



Fig. 8. — CO(2— 1) contours at intermediate resolution overlaid on a selection of multiwavelength data. The labels in the top left corners 
refer to the greyscale images - see text for details. 

clumps are entirely composed of molecular gas and using 
-Mviriai/-^co(2-i) as an es ti ma te of aco, we find aco for 
the clumps is <0.2-0.7 M (K km s" 1 pc 2 )" 1 . 

Recently a couple of other observational studies have 
looked at molecular gas clumps in z > 1 sta r forming 
galaxies. For example, iTacconi et al.l (|2010[ ) observed 
molecular gas clumps in a 'normal' star forming galaxy at 
z <~ 1.2. These clumps had gas masses of — 5xl0 9 Mq, 
diameters <2-4 kpc, brightness temperatures of ~10- 
25 K, and gas volume densities ~100 cm -3 , similar to 
what we get for GN20 if we assume that the clumps are 
roughly the beamsize. Their velocity dispersions, how- 
ever, were only ~20 km s , causing them to postulate 
that the 'clumps' observed were more likely loose con- 
glomerates of several GMCs. Another recent paper by 
ISwinbank et al.l (|2011[) identified molecular gas clumps in 
a lensed star forming galaxy at z = 2.3, SMM J2135. The 
clumps in SMM J2135 had gas masses of -0.3-1 xlO 9 
Mq, a median diameter of 200 pc, and gas volume densi- 
ties of ~1, 000-40, 000 cm~ 3 . To achieve volume densities 
this high, the clumps in GN20 would need to be up to an 
order of magnitude smaller than the current beamsize. 

4.3. Multiwavelength Comparison 

Figure [5] shows the 0th moment CO (2-1) contours 
at intermediate resolution overlaid on a selection of 
the multiwavelength data available fo r GN20, includ- 



where C = 5 is the dimensionless prefactor for a uniform 
sphere, a v is the observed ID velocity dispersion in km 
s _1 , R g is the graviational radius, and G — 1/232 pc 
(km/s) 2 Mq 1 is the gravitational constant. If we take 
the beam HWHM as an upper limit on the gravitational 
radius, we estimate virial masses for the clumps of < 1.0- 
3.6 x 10 9 M Q (Table [2]). These masses are all within 
l-2cr of the observed molecular gas masses, even assum- 
ing a low ULIRG-like mass conversion factor of aco = 
0.8 Mq (K km s _1 pc 2 ) -1 and ignoring the possible pres- 
ence of significant stellar mass. Within the considerable 
uncertainties, it is possible that the clumps have masses 
consistent with self-gravitation. 

The gas masses derived for the clumps are already on 
the high side (in general) of the virial mass estimates, 
even with a low conversion factor and despite the fact 
that the virial mass estimates are upper limits. This 
implies that for the majority of the clumps, the gas 
masses would significantly overshoot the virial masses 
for aco >> 0.8 M Q (K km s" 1 pc 2 )^ 1 . Assuming the 



ing th e HST+ACS z850-band (top left: iGiavalisco et al.l 
12001 . the WFC3 140W-band (top right), the WIRCAM 
K-band (bottom left), and the IRAC 3.6/im data 
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Fig. 9. — Spatially-resolved molecular gas excitation in GN20. The left panels show the CO(2-l) and CO(6-5) images at the same 
resolution (0.84" X 0.67") and beam position angle. The right plot shows the result of dividing the CO(6— 5) image by the CO(2— 1) image 
(error bars do not include overall flux calibration uncertainties). Noise dominates the division at radii larger than ~1.2". The dotted and 
dashed lines are models (at the same resolution) where the CO(6-5) emission profile has an exponential scale length that is 2 and 5 times 
more compact than the CO (2—1), respectively. These model images have been normalized by the peak value in the image and are meant 
to demonstrate that changes in the excitation ratio are detectable at this resolution. 

(bottom right). Not shown is the high er-resolution 
VLA+MERLIN data (ICasev et all 120091) whic h was 
compared to GN20 in CO in ICarilli et all (|2010l ). The 
object detected in the HST+ACS z850-band data (top 
left panel) was followed up spectroscopically with the 
Deep Imaging Multi-Object Spectrograph (DEIMOS) on 
Keck II and determ ined to have a redshift of z=4.055 
(jDaddi et a l. 2009b), consistent with the redshift derived 
for GN20 from its CO lines. 

Assuming the object is related to GN20, as indicated 
by the good match in redshift, then the question is 
what is causing the clear offset from the CO and ra- 
dio/ (sub) millimeter posi tion also reported in previous 
work on GN20 (for e.g.. ICarilli et al.l[20Tot iPope et all 
[20061: lYounger et all 120081: llono et al.M2006l) . We have 
confirmed that the astrometry of the HST imaging is 
good to <0.15" (IDaddi et alll2007fl . so astrometric error 
is unlikely. In addition, the counterpart in the WFC3 
140W-band imaging is similar in position/morphology. 
Such large offsets between dust /molecular gas emission 



and (rest-frame) U V/ optical imaging have been seen in 
a number of SMGs dChapman et al l 120051 : iCapak et all 



120081: iRiechers et al.N2010D . A possible hint in this case 
comes from the WIRCAM K-band image, which is sim- 
ilarly offset, but shows extended emission in the direc- 
tion of the radio/ (sub) millimeter counterparts. This ex- 
tended emission is significant at the ~6cr level. A pro- 
gression through the observations shows that the peak 
of emission shifts toward the radio position as we sam- 
ple longer wavelengths. One interpretation is therefore 
that the radio/ (sub)millimeter emission is offset from 
the (observed) optical emission because the dust that 
is presumably associated with the molecular gas heav- 
ily obscures the rest-frame UV and optical emission. A 
significant amount of dust obscuration is not surprising 
in a starburst of several thousand M yr" 1 , the SFR 



estima ted for GN20 (|Magdis et all 120111 : IDaddi et all 
l2009bf ). Heavy dust obscuration is thought to explain 
the complet e absence of optica l counterparts for some 
SMGs (e.g.. iWalter et all 120121 ) . The top left panel in 
Figure [S] shows that the rest-frame UV emission is well 
aligned with the edge of the 0th moment CO contours, 
giving the clearest indication that the ACS counterpart is 
physically related but is simply offset due to differential 
dust obscuration. 

4.4. Spatially-resolved Gas Excitation 

To investigate the spatially-resolved molecular gas ex- 
citation, we compare our CO(2-l) data with the CO(6- 
5) Plateau d e Bu re Interferometer (PdBI) data from 
ICarilli et all (|2010f ). The CO(6-5) data are uniformly- 
weighted at an angular resolution of 0.84" x 0.67". We 
tapered the CO(2-l) data to the same resolution and 
imaged over the same velocity range (800 km s _1 ). 

We then used the AIPS task IRING to calculate the 
average CO(6-5)/CO(2-l) ratio as a function of radius. 
The centroids of the CO(6-5) and CO(2-l) emission 
peaks were consistent within 0.08", so we used the aver- 
age of the positions as the center of the concentric annuli. 
We set the width of the annuli to be 0.25", i.e., one-third 
of the synthesized beam, and we specified a major axis 
position angle of 115° and an inclination of 30° as deter- 
mined in Section 14.1.21 

The results are shown in Figure |H1 The left panels 
show the CO(2-l) and CO(6-5) emission, using the na- 
tive beamshape of the CO (6-5) map. The plot on the 
right is the result of dividing the CO(6-5) emission by 
the CO(2-l) emission, averaging in elliptical annuli as 
explained above. Note that the error bars are purely 
statistical and do not account for the overall flux cali- 
bration uncertainties in the maps (which in both cases 
are of order 10%). 
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The CO(6-5)/CO(2-l) ratio is an indicator of the exci- 
tation state of the gas. The ratio hovers around a value 
of —3, a factor of three lower than expected for ther- 
mal excitation. This indicates that the CO(6-5) line is 
sub-thermally e xcited, consis t ent w ith the CO excitation 
ladder shown in lCarilli et al.l (|2010t) . Carilli ct al. deter- 
mined that the data were best fit with a two-component 
gas model including a diffuse, low-excitation component 
and a more concentrated high-excitation component. 

With the data at hand, we can go one step further and 
investigate the radial dependence of the CO ratio. From 
Figure [9l we see that the ratio is consistent with a con- 
stant value out to at least 1", corresponding to the ob- 
served radius of the CO(2-l) emission in the Oth moment 
map. Note that we have spaced the radial bins so as to 
achieve Nyquist sampling, so the first two data points are 
within the beam. Nevertheless, this result tells us that 
to first order, the CO(6-5) emission has the same radial 
profile as the CO(2-l) emission. This would not be the 
case if, for example, the CO (6-5) emission (which is more 
closely related to star formation) was more tightly con- 
centrated than the lower-excitation CO (2-1) emission. 
To illustrate this point, we show two different models 
(at the same resolution as the data) where the CO(6-5) 
emission is more centrally-concentrated than the CO (2- 
1) emission. In the first model (dotted line), the expo- 
nential scale length of the CO(6-5) emission is two times 
smaller than that of the CO (2-1) emission. In the sec- 
ond model (dashed line), the CO(6-5) emission has a 
scale length that is five times smaller. In both cases, 
the CO (2-1) disk parameters come from our dynamic 
modeling in Section 14.1.21 The models have been nor- 
malized to the peak value in the maps and are meant 
to demonstrate the clear drop-off in the excitation ratio 
that would be detectable for a central star forming event 
even at this resolution. Instead, the flatness we observe 
is consistent with the idea that the two transitions are 
cospatial, and that the star formation is therefore spread 
out over a large portion of the (14 kpc diameter) disk. 

5. WHAT IS DRIVING THE STARBURST? 

GN20 has several properties which make it unique as 
an SMC It has a higher redshift than the typical SMC, 
and it is extremely bright - its 850 fim flux density of 20.3 
mJy m akes it one of the m ost luminous starburst galaxies 
known (jPope et al.l 120061 ). In addition, it appears to be 
part of a massive proto-cluster of galaxies: there are two 
other z ~ 4.05 SMGs ju st — 25"/170 kpc away known as 
GN20.2a and GN20.2b (IDaddi et al.ll2009bt [Carilli et al.1 
IMll : there are 14 z ~ 4 Lyman Break Galaxies within 
25"; and there is an overdensity of z > 3.5 IRAC selected 
galaxies in the field. 

The question is: what process is driving the gas accre- 
tion and conversion into stars in this example of one of 
the earliest extreme starburst galaxies? The majority of 
SMGs are thought to be gas-rich major mergers. There 
are multiple examples of SMGs with very disturbed kine- 
matics and even showing multiple components, provid- 
ing d irect observational evidence for the merger picture 
(e"g. iTacconi et al.ll2008UEngel et al.ll2010t Hvison et al.l 
[20lil: iRiechers et alJl2011aUbD For example, by identify- 
ing likely radio and/or mid-IR ide ntifications to SMG s 
in the SHADES S ource Catalog (ICoppin et al.l l2006h . 
Ilvison et~a l. (2007) found that significantly more SMGs 



have multiple robust counterparts than would be ex- 
pected by chance alone. 

In GN20, on the other hand, we observe an extended, 
rotating gas disk with a velocity field that is consistent 
with a flat rotation curve. The gas reservoir appears 
to be fragmented into multiple clumps, with each clump 
making up only a few percent of the total mass. This 
morphology stands in contrast to prototypical merging 
sources like Arp 220 and other ULIRGs, where the molec- 
ular gas is typically concentrated in one or two dense, 
centrally-located regions. If the diffuse and clumpy 
gas components correspond to the two gas excitation 
states predicted to exist in GN20 with LVG modeling 
(|Carilli et al.l [2010^). then the clumpy component (pre- 
sumably given by the higher excitation state gas) only 
makes up —half of the total gas mass, with the remain- 
der existing as extended, diffuse emission. 

Recently, cold mode accretion (CMA) has been sug- 
gested as a possible alternative to the major merger 
formation scenario for some SMGs flFardal et all [2 001: 
iFinlator et al.ll2006tlDave et al.l2O10HCarilli et alJ2 010). 
In this scenario, the star formation in these SMGs 
would be instead fueled by minor mergers and the 
smooth infall of cold gas along streams. Unlike a 
major merger, this process can leave the disk rela- 
tively undi sturbed, especiall y for a galaxy as massive 
as GN20 (|Dave et al.l I2010D . However, wet mergers 
have been shown to rapidly relax into smooth disks, 
so the presence of ordered rotation alone cannot im- 
mediately be taken as an argument against a ma- 
jor merger dRobertson et alj|2006t iRobertson fc Bullock! 
12005 iBouTnaud fc Elmegreenll2009ft .~ 

If GN20 were a major merger where the gas had re- 
virialized into a rotating disk, one may ask whether 
it would still be an SMG at this late stage? Accord- 
ing to hydrodynamic simulations of major mergers, the 
submm-bright phase in a (z — 2) major merger-driven 
starburst should occur as the gal axies ap proach final coa- 
lescence dNaravanan et al.ll2009t but c.f. lHavward et al.l 
l20ll 120121). During this brief phase (-0.03 Gyr), the 
850 (im flux is significantly boosted and the galaxy may 
be selected as a luminous submillimeter source. How- 
ever, according to the simulations, this phase coincides 
with the stage when the disk-like morphology is most 
disturbed, with the large st fraction of gas tida lly dis- 
rupted from disk rotation (|Naravanan et a l. 2009, Figure 
6). These simulations also show that the gas during the 
SMG phase is relatively concentrated, with a character- 
istic radius of —1.5 kpc. Such a compact morphology 
is expected, since the large-scale gravitational torques 
induced by gas-rich major mergers are efficient at re - 
moving angular momentum (B arnes fc Hernquisdll996T ). 
In GN20, on the other hand, we simultaneously observe 
1) an extremely bright (Ssso = 20.3 mJy) submillimeter 
flux, 2) ordered rotation, and 3) a very extended (14 kpc) 
gas disk. Taken together, these observations may be an 
indication that the star formation in GN20 is fueled by 
some process other than a major merger. A simi lar argu- 
ment was put forward by iBothwell et al.l (|2010| ) regard- 
ing the SMG HDF132. Of course, this argument rests 
on the assumption that these simulations provide an ac- 
curate representation of merger-driven SMG formation 
for all SMGs, which may not be the case, especially for 
notable SMGs like GN20. 
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On the other hand, there are certain properties of 
GN20 which may be difficult to explain without a ma- 
jor merger. In particular, the rotational velocity derived 
from the dynamical modeling is extremely high (576 km 
s _1 ). While we cautioned that this is largely dependent 
on the galaxy's inclination, (which is hard to constrain at 
our S/N), a rotational velocity this large may be hard to 
reproduce in a spinning disk embedded in its dark matter 
halo, and the linewidth could be more easily explained 
as resulting from the final coalescence stage of a major 
merger. Also noteworthy is the observation that GN20 is 
undergoing an extremely dusty mass ive starburst. This i s 
supported by its high specific S FR (iDaddi et al.ll2"0 09b N ). 
and also by its hot FIR SED (|Magdis et all 120111 ) and 
low CO-to-H2 conversion factor, which both imply a 
high star formation efficiency. The critical question is: 
is it possible to form such an object without relying on a 
major merger? Locally, at least, the answer is no - such 
objects are all major mergers (ULIRGs). 

Perhaps the real answer lies s o mewh ere in between 
the two scenarios. iFinlator et a l. (2006) studied z ~ 4 
B-band dropout galaxies in cosmological hydrodynamic 
simulations and found two instances of galaxies form- 
ing stars at over 1000 M Q yr _1 , several times that of 
their 100 Myr averages. These objects both resided in 
highly overdense regions resembling the cores of massive 
clusters. The authors found that, while a merger of two 
massive galaxies was not the dominant mechanism for fu- 
eling star formation, there were still signs of interaction 
with many of the smaller nearby galaxies. GN20 lies in a 
similarly dense environment - the overdensity of B-band 
dropouts and IRAC-selected galaxies suggests a proto- 
cluster of total mass 10 14 M© centered o n GN20 and 
the ne arby SMGs GN20.2a and GN20.2b (jDaddi et al.l 
2009b). It is likely, therefore, that even if there is a dom- 
inant rotating disk, enhanced star formation over the disk 
is being triggered by the gravitational torques from the 
other galaxies. 

6. CONCLUSIONS 

We have presented VLA observations of CO(2-l) in 
the z = 4.05 SMG GN20. These high-resolution data 
allow us to do detailed imaging of the molecular gas on 
scales down to 1.3 kpc. We use these data to model 
the overall gas dynamics and determine the properties of 
individual star forming clumps in an extreme starburst 
galaxy just 1.6 Gyr after the Big Bang. We summarize 
our main findings in the following. 

The data reveal a clumpy, extended gas reservoir in 
unprecedented detail. A carefully constructed 0th mo- 
ment map recovers emission over an area 14 ± 4 kpc in 
diameter, and the first moment map shows a clear veloc- 
ity gradient. A dynamical analysis shows that the data 
are consistent with a large, rotating disk with a maxi- 
mum rotational velocity of v max — 575 ± 100 km s _1 , 
an inclincation of 30° ± 15°, and a dispersion of 100 ± 
30 km s _1 . While it is difficult to determine the exact 
shape of the rotation curve, a flat rotation curve gives 
the best fit to the data. 

From the dynamical modeling, we derive a dynamical 
mass for GN20 of 5.4 ± 2.4 x 10 11 M . We use this 
dynamical mass in combination with the stellar mass to 
put constraints on the CO-to-Bi2 mass conversion factor 
(a co ), finding a C o = 1-1 ±0.6 M (K km s" 1 pc 2 )" 1 . 



This value is consistent with a low, ULIRG-like value, 
although better data is needed to confirm this result. Our 
estimate is also in agreement with recent limits set for 
GN20 using the local M gas /Ma U st vs. metallicity relation. 

We use two different source-finding techniques to iden- 
tify five distinct molecular gas clump candidates on the 
disk of GN20. These clumps have masses of a few percent 
of the total gas mass and a median linewidth of 140 km 
s" 1 (FWHM). Their surface densities are >3, 200-4,500 
x (cfCo/0.8) Mq pc -2 , corresponding to volume densi- 
ties of >100 cm~ 3 (where the lower limit corresponds to 
GMCs in the Milky Way). Their peak brightness temper- 
atures (16-31 K) are approaching both the average dust 
temperature (33K) and the 'warm' component of a two- 
component LVG model, implying that we may be close to 
resolving the clumps. Within the substantial uncertain- 
ties, all clumps are consistent with being self-gravitating 
masses, and we constrain their CO-to-Bi2 mass conver- 
sion factor to be <0.2-0.7 M (K km s _1 pc 2 )" 1 . 

We compare the observed distribution of CO to the 
emission seen at other wavelengths. This multiwave- 
length comparison demonstrates that the molecular gas 
is concentrated in a region of the galaxy which is heavily 
obscured in the rest-frame UV and optical. The optical 
counterpart likely constitutes a small percentage of the 
stellar light. 

Using previous observations of the CO(6-5) in GN20, 
we examine the CO(6-5)/CO(2-l) excitation ratio in 
elliptically-averaged annuli. The excitation ratio is con- 
sistent with a constant value over the extent of the 14 
kpc disk. This result implies that the star formation is 
spread out over a large portion of the disk, rather than 
concentrated in a central star forming event. 

We discuss our results in the context of different fueling 
scenarios for GN20. While the presence of an extended, 
rotating disk during the SMG phase may point toward a 
process other than a major merger (e.g., cold mode accre- 
tion), its large linewidth and extremely dusty starburst 
would be more easily explained in a merger picture. If 
the star formation is fueled primarily by CMA, then, due 
to GN20's location in a dense protocluster environment, 
it is likely that it is being enhanced by interactions with 
nearby galaxies. 

We conclude by highlighting the importance of low-J 
CO studies, which trace the bulk of the molecular gas, 
in studying the formation of massive galaxies throughout 
cosmic time. Such studies cannot currently be done in a 
statistical way due to the time required to attain suffi- 
cient S/N at high-resolution, even with state-of-the-art 
instruments such as the VLA. This fact makes targeted 
case studies such as this that much more important. 
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